package scu.maqiang.homogenization;

import scu.maqiang.fes.*;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class SOTS3DT4102 {
	public SOTS3DT4102(FES3T42 fsT4, FES3T102 fsT10) {
		this.fsT4 = fsT4;
		this.fsT6 = fsT10;
	}

	public void elasticityCellSolution(MatrixFunc conFunc, double[][][] param, BVPType tp, double[][] FOCS,
			double[][] D, double[][] SOCS) {
		SRMatrix K = new SRMatrix(fsT4.GetNdof());
		fsT4.assembleStiff(conFunc, param, tp, K);
		double[][] RHS = new double[6][fsT4.GetNdof()];
		fsT4.assembleGlobalVector(conFunc, param, SOTS3DT4102::elementElasticityFirstOrder, tp, RHS);
		fsT4.applyBC_MBN(K, Direct.All, 1, 2, 3, 4, 5, 6);
		NewIterSSolver solver = new NewIterSSolver(K);
		MVO.fill(FOCS, 0.0);
		for(int i = 0; i < 6; i++) {
			fsT4.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
			solver.PCGSSOR(RHS[i], FOCS[i], 1.75, 1); // 11
		}

		computeHomogenizedElasticity(conFunc, param, FOCS, D);
		
		//System.out.println(MVO.toString(D));
		//System.exit(0);
		RHS = new double[18][fsT4.GetNdof()];
		fsT4.assembleGlobalVector(conFunc, param, FOCS, D, SOTS3DT4102::elementElasticitySecondOrder, tp, RHS);
		MVO.fill(SOCS, 0.0);
		for(int i = 0; i < 18; i++) {
			fsT4.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
			solver.PCGSSOR(RHS[i], SOCS[i], 1.75, 1);
		}
	}

	public static void elementElasticityFirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleXY) {
		double[][] grad = new double[4][3];
		double vol = CG.gradient3DT4(xy, grad);
		
		double[][] B = new double[12][6];
		B[0][0] = B[1][3] = B[2][5] = -grad[0][0] * vol;
		B[1][1] = B[0][3] = B[2][4] = -grad[0][1] * vol;
		B[2][2] = B[1][4] = B[0][5] = -grad[0][2] * vol;

		B[3][0] = B[4][3] = B[5][5] = -grad[1][0] * vol;
		B[4][1] = B[3][3] = B[5][4] = -grad[1][1] * vol;
		B[5][2] = B[4][4] = B[3][5] = -grad[1][2] * vol;

		B[6][0] = B[7][3] = B[8][5] = -grad[2][0] * vol;
		B[7][1] = B[6][3] = B[8][4] = -grad[2][1] * vol;
		B[8][2] = B[7][4] = B[6][5] = -grad[2][2] * vol;

		B[9][0] = B[10][3] = B[11][5] = -grad[3][0] * vol;
		B[10][1] = B[9][3] = B[11][4] = -grad[3][1] * vol;
		B[11][2] = B[10][4] = B[9][5] = -grad[3][2] * vol;

		eleXY[0] = MVO.matmul(B, coef[0]);
		eleXY[1] = MVO.matmul(B, coef[1]);
		eleXY[2] = MVO.matmul(B, coef[2]);
		eleXY[3] = MVO.matmul(B, coef[3]);
		eleXY[4] = MVO.matmul(B, coef[4]);
		eleXY[5] = MVO.matmul(B, coef[5]);
	}

	public void computeHomogenizedElasticity(MatrixFunc conFunc, double[][][] param, double[][] FOCS, double[][] D) {
		double[][] grad = new double[4][3];
		double[][] xy = new double[4][5];
		Mesh3T4 mesh = (Mesh3T4) fsT4.getMesh();
		double[][] N = new double[6][12];
		double[][] B = new double[6][12];
		double wholeVol = 0.0;
		MVO.fill(D, 0.0);
		for (int k = 0; k < mesh.getNt(); k++) {
			int[] ele = mesh.getElement(k);
			mesh.getCoordPerElement(ele, xy);
			int[] dof = fsT4.getDof(ele);
			MVO.getSlice(FOCS[0], dof, N[0]);
			MVO.getSlice(FOCS[1], dof, N[1]);
			MVO.getSlice(FOCS[2], dof, N[2]);
			MVO.getSlice(FOCS[3], dof, N[3]);
			MVO.getSlice(FOCS[4], dof, N[4]);
			MVO.getSlice(FOCS[5], dof, N[5]);
			double vol = CG.gradient3DT4(xy, grad);

			B[0][0] = B[3][1] = B[5][2] = grad[0][0];
			B[1][1] = B[3][0] = B[4][2] = grad[0][1];
			B[2][2] = B[4][1] = B[5][0] = grad[0][2];

			B[0][3] = B[3][4] = B[5][5] = grad[1][0];
			B[1][4] = B[3][3] = B[4][5] = grad[1][1];
			B[2][5] = B[4][4] = B[5][3] = grad[1][2];

			B[0][6] = B[3][7] = B[5][8] = grad[2][0];
			B[1][7] = B[3][6] = B[4][8] = grad[2][1];
			B[2][8] = B[4][7] = B[5][6] = grad[2][2];

			B[0][9] = B[3][10] = B[5][11] = grad[3][0];
			B[1][10] = B[3][9] = B[4][11] = grad[3][1];
			B[2][11] = B[4][10] = B[5][9] = grad[3][2];
			double[] ip = mesh.getinteriorPointInElement(xy);
			double[][] con = conFunc.action(ip, mesh.getElementLabel(k), param);
			for (int i = 0; i < 6; i++) {
				for (int j = 0; j < 6; j++) {
					D[i][j] += (con[i][j] + MVO.dot_product(con[i], MVO.matmul(B, N[j]))) * vol;
				}
			}
			//System.out.println("OK " + k + "of" + mesh.getNt());
			wholeVol += vol;
		}
		MVO.mul(D, 1.0 / wholeVol);
	}

	public static void elementElasticitySecondOrder(int[] ele, double[][] xyz, double[][] con, double[][] FOCS, double[][] D,
			BVPType tp, double[][] eleXY) {
		double[][] grad = new double[4][3];
		double vol = CG.gradient3DT4(xyz, grad);
		double vol4 = 0.25 * vol;
		double[][] B = new double[6][12];
		B[0][0] = B[3][1] = B[5][2] = grad[0][0];
		B[1][1] = B[3][0] = B[4][2] = grad[0][1];
		B[2][2] = B[4][1] = B[5][0] = grad[0][2];

		B[0][3] = B[3][4] = B[5][5] = grad[1][0];
		B[1][4] = B[3][3] = B[4][5] = grad[1][1];
		B[2][5] = B[4][4] = B[5][3] = grad[1][2];

		B[0][6] = B[3][7] = B[5][8] = grad[2][0];
		B[1][7] = B[3][6] = B[4][8] = grad[2][1];
		B[2][8] = B[4][7] = B[5][6] = grad[2][2];

		B[0][9] = B[3][10] = B[5][11] = grad[3][0];
		B[1][10] = B[3][9] = B[4][11] = grad[3][1];
		B[2][11] = B[4][10] = B[5][9] = grad[3][2];
		double[][] BT = MVO.transpose(B);

		double[][] w = { { vol4, 0.0, 0.0}, { 0.0, vol4, 0.0}, {0.0, 0.0, vol4}, { vol4, 0.0, 0.0}, { 0.0, vol4, 0.0}, {0.0, 0.0, vol4}, 
				{ vol4, 0.0, 0.0}, { 0.0, vol4, 0.0}, {0.0, 0.0, vol4}, {vol4, 0.0, 0.0}, {0.0, vol4, 0.0}, {0.0, 0.0, vol4}}; //
		double[][] wT = MVO.transpose(w);
//
		double[] C = new double[3];
		double[] part1 = null;
		double[][] con2 = new double[3][];
		double[][] con2T = null;
		double[] part2 = null;
		
		for(int k = 0; k < 6; k++) {
			
			C[0] = con[0][k] - D[0][k] + MVO.dot_product(con[0], MVO.matmul(B, FOCS[k]));
			C[1] = con[3][k] - D[3][k] + MVO.dot_product(con[3], MVO.matmul(B, FOCS[k]));
			C[2] = con[5][k] - D[5][k] + MVO.dot_product(con[5], MVO.matmul(B, FOCS[k]));
			part1 = MVO.matmul(w, C);

			con2[0] = con[0];
			con2[1] = con[3];
			con2[2] = con[5];
			con2T = MVO.transpose(con2);
			part2 = MVO.matmul(MVO.matmul(MVO.matmul(BT, con2T), wT), FOCS[k]);
			
			//System.out.println(MVO.toString(part1));
			//System.out.println(MVO.toString(part2));
			
			for (int i = 0; i < 12; i++) {
				eleXY[3 * k][i] = part1[i] - part2[i]; // 111
			}
			
			//System.out.println(MVO.toString(eleXY[0]));
			//System.exit(0);
			
			//=======================================================
			C[0] = con[3][k] - D[3][k] + MVO.dot_product(con[3], MVO.matmul(B, FOCS[k]));
			C[1] = con[1][k] - D[1][k] + MVO.dot_product(con[1], MVO.matmul(B, FOCS[k]));
			C[2] = con[4][k] - D[4][k] + MVO.dot_product(con[4], MVO.matmul(B, FOCS[k]));
			part1 = MVO.matmul(w, C);

			con2[0] = con[3];
			con2[1] = con[1];
			con2[2] = con[4];
			con2T = MVO.transpose(con2);
			part2 = MVO.matmul(MVO.matmul(MVO.matmul(BT, con2T), wT), FOCS[k]);
			for (int i = 0; i < 12; i++) {
				eleXY[3 * k + 1][i] = part1[i] - part2[i]; // 121
			}
			
			//=======================================================
			C[0] = con[5][k] - D[5][k] + MVO.dot_product(con[5], MVO.matmul(B, FOCS[k]));
			C[1] = con[4][k] - D[4][k] + MVO.dot_product(con[4], MVO.matmul(B, FOCS[k]));
			C[2] = con[2][k] - D[2][k] + MVO.dot_product(con[2], MVO.matmul(B, FOCS[k]));
			part1 = MVO.matmul(w, C);

			con2[0] = con[5];
			con2[1] = con[4];
			con2[2] = con[2];
			con2T = MVO.transpose(con2);
			part2 = MVO.matmul(MVO.matmul(MVO.matmul(BT, con2T), wT), FOCS[k]);
			for (int i = 0; i < 12; i++) {
				eleXY[3 * k + 2][i] = part1[i] - part2[i]; // 131
			}
		}
		
		//System.out.println(MVO.toString(eleXY));
		//System.exit(0);

	}

	public void assembleElasticitySOTSSolution(FES2T32 fsFine, double[] T_0, double[] T_1, double[] T_2,
											   FES2T62 fsHomo, double[] T_Homo, FES2T32 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		
	}

	FES3T42 fsT4;
	FES3T102 fsT6;
}
